Supplement A: Estimating Travel Cost

Published

January 23, 2024

1 Overview

Goal: estimate cost-distance (in time) from each ranch to Ciudad Constitucion, La Paz, and local springs.

Data: a Digital Elevation Model (DEM) along with point locations for ranche clusters, springs, cities.

Method: Campbell’s hiking function (Campbell et al. 2019).

There are two watersheds in the project area. Each one has a single road that connects it to Highway 1 and thus to the markets in Ciudad Constitucion to the north and the capitol La Paz to the south. The “road” that connects the watersheds to each other is not maintained and thus rarely used. As a consequence, nearly all of the variance in travel time to cities is found within each watershed. For that reason, we estimate the time it takes to leave each watershed (“leaving” here means hitting the edge of the project area), rather than the true time it takes to get to each city. We interpret this as the relative difference in time it takes to get to Ciudad Constitucion from the northern watershed and the relative difference in time it takes to get to La Paz from the southern watershed. To get to La Paz from the northern watershed, we simply add an arbitrary half hour to the total time to get out of that watershed. The same would presumably be the case to get to Ciudad Constitucion from the southern watershed, though that value does not factor into our analysis.

For springs, Dr. MacFarlan’s ethnographic data suggest that individuals residing within the same ranch cluster tend to rely on one or two springs at most. Thus, we estimate road distance to the two nearest springs and take the average.

All spatial data are projected to the Mexico ITRF92 / UTM zone 12N CRS (EPSG:4485).

2 R Preamble

library(dplyr)
library(ggnewscale)
library(ggplot2)
library(here)
library(igraph)
library(sf)
library(terra)
library(viridis)

Specify consistent plot theme here:

Code
# this is designed to position the legend 
# at the top left corner of the plot area

theme_set(theme_bw(12))

theme_update(
  axis.text = element_text(size = rel(0.5), color = "gray45"),
  axis.text.y = element_text(angle = 90, hjust = 0.5),
  axis.ticks = element_line(color = "gray45", linewidth = 0.1),
  axis.ticks.length = grid::unit(0.1, "cm"),
  axis.title = element_blank(),
  legend.background = element_rect(fill = "transparent", color = "transparent"),
  legend.key.height = grid::unit(0.4, "cm"),
  legend.key.width = grid::unit(0.3, "cm"),
  legend.margin = margin(l = 0),
  legend.justification = c("left", "top"),
  legend.spacing.x = grid::unit(0.05, "cm"),
  legend.spacing.y = grid::unit(0.05, "cm"),
  legend.text = element_text(vjust = 0.7, margin = margin()),
  legend.title = element_blank(),
  panel.grid = element_blank(),
  plot.title = element_text(size = rel(1.1), face = "bold")
)

point_colors <- c("#FDE74C", "#EE7733", "#0077BB")

pretty_breaks <- function(x, .axis = "x", dx = 2500, n = 4){ 
  
  mn <- paste0(.axis, "min")
  mx <- paste0(.axis, "max")
  
  brks <- seq(x[[mn]] + dx, x[[mx]] - dx, length.out = n)
  
  round(brks, digits = -3)
  
}

3 Data

3.1 Geopackage

Specify path to geopackage database holding all spatial vector data.

gpkg <- here("data", "choyero.gpkg")

3.2 Project Window

bcs <- read_sf(gpkg, layer = "baja")

window <- read_sf(gpkg, layer = "window")

roads <- read_sf(gpkg, layer = "roads")

watersheds <- read_sf(gpkg, layer = "watersheds")

3.3 Elevation (DEM)

dem <- rast(here("data", "dem_30m.tif"))

3.4 O-D Points

# ORIGIN ---
clusters <- gpkg |> 
  read_sf(layer = "clusters") |>
  mutate(variable = "cluster")

# DESTINATION --- to city, terminus at edge of window
terminus <- gpkg |> 
  read_sf(layer = "terminus") |>
  mutate(variable = "terminus")

springs <- gpkg |> 
  read_sf(layer = "springs") |>
  select(id, name) |>
  mutate(variable = "spring")
Code
od_points <- clusters |> 
  rename("name" = cluster) |> 
  mutate(name = as.character(name)) |> 
  select(name, variable) |> 
  bind_rows(terminus, springs) |> 
  mutate(variable = ifelse(variable == "terminus", "city", variable))

bb8 <- st_bbox(window)
dx <- bb8[["xmax"]] - bb8[["xmin"]]
dy <- bb8[["ymax"]] - bb8[["ymin"]]

# add hillshade relief to visualize elevation
# https://dominicroye.github.io/en/2022/hillshade-effects/
hillshade <- shade(
  slope = terrain(dem, "slope", unit = "radians"),
  aspect = terrain(dem, "aspect", unit = "radians"),
  angle = 45,
  direction = seq(45, 360, by = 45),
  normalize = TRUE
)

hillshade <- mask(sum(hillshade), vect(st_union(watersheds)))
names(hillshade) <- "shade"

gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hillshade, xy = TRUE),
    aes(x, y, fill = shade)
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide = "none"
  ) +
  new_scale_fill() +
  geom_raster(
    data = as.data.frame(dem, xy = TRUE),
    aes(x, y, fill = elevation),
    alpha = 0.7
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide="none"
  ) +
  new_scale_fill() +
  geom_sf(
    data = watersheds, 
    fill = "transparent", 
    color = alpha("darkblue", 0.35),
    linewidth = 0.2
  ) +
  geom_sf(
    data = roads, 
    color = "black",
    linewidth = 0.3
  ) +
  geom_sf(
    data = od_points, 
    aes(shape = variable, fill = variable, color = variable),
    stroke = 0.3,
    size = 1.8
  ) +
  scale_color_manual(values = c("grey15", "white", "white")) +
  scale_fill_manual(values = alpha(point_colors, 0.75)) +
  scale_shape_manual(values = c(23, 22, 21)) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y", n = 3), 
    expand = expansion(add = 1000)
  ) +
  theme(legend.position = c(0.03, 0.98))

ggsave(
  plot = gg,
  filename = here("figures", "od-points.png"),
  dpi = 800,
  width = 3.5,
  height = 3.5 * (dy/dx)
)

4 Path analysis

To estimate least-cost paths, we first derive slope estimates from the DEM then calculate the hiking speed at which individuals can traverse any given cell using Campbell’s hiking function. After that, we decrease the cost in units of time to traverse grid cells crossed by roads and other paths using weights that approximate travel times by car or horse. This is equivalent to increasing the speed one travels in those cells. Below are the specific weights we use.

Down-weight for travel-cost
Type Weight Approx. Speed
primary 0.2 25 mph (41.5 kmh)
secondary 0.5 10 mph (16.2 kmh)
tertiary 0.8 07 mph (11.0 kmh)

The necessary functions are defined in the folded code chunk below.

Code
#' Cost surface
#'
#' Create a `terrain` representing the cost of travel. 
#'
#' The `SpatRaster` must have a projected CRS, with both distance units and elevation
#' values in meters.
#'
#' @param x a `SpatRaster` with elevation values.
#' @param hf a character string specifying the preferred hiking function. Options
#'   include `tobler` and `campbell` (default).
#' @param max_slope a numeric value specifying the maximum allowed hiking slope
#'   _in degrees_ (default is 45). Suppose you ask: Would a person really hike
#'   _that_? This is the shallowest slope where the answer to that question is No.
#' @param neighbors numeric, one of 4, 8, or 16 (default) specifying the number of
#'   adjacent cells between which to calculate movement cost. These are passed on
#'   to `terra::adjacent()` and represent standard neighborhoods (rook, queen, and
#'   knight + queen).
#' @param ... Arguments passed on to the hiking function specified by `hf`.
#'
#' @return a `terrain` representing cost of travel.
hf_terrain <- function(x,
                       hf = "campbell",
                       max_slope = 45,
                       neighbors = 8,
                       ...) {

  # find adjacent cells in the raster
  # make sure to ignore all na values
  n_cells <- terra::ncell(x)
  
  na_cells <- which(is.na(terra::values(x)))
  
  cells <- setdiff(1:n_cells, na_cells)
  
  adj <- terra::adjacent(
    x,
    cells = cells,
    directions = neighbors,
    pairs = TRUE
  )
  
  adj <- adj[!adj[, 2] %in% na_cells, ]
  
  ### slope ###
  # change in y between adjacent pixels
  heights <- terra::values(x)[, 1]
  
  rise <- (heights[adj[, 2]] - heights[adj[, 1]])
  
  # change in x between adjacent pixels
  run <- terra::distance(
    terra::xyFromCell(x, adj[, 1]),
    terra::xyFromCell(x, adj[, 2]),
    lonlat = FALSE,
    pairwise = TRUE
  )
  
  slope <- rise/run
  
  ### speed ###
  speed_fun <- switch(
    hf,
    "campbell" = campbell,
    "tobler"   = tobler,
    stop(paste0("hiker does not implement the ",
                hf,
                " hiking function."))
  )
  
  speed <- speed_fun(slope, ...)
  
  ### conductance ###
  conductance <- speed/run
  
  # what is a realistic slope people would try to hike?
  max_slope <- tan(max_slope * pi / 180) # degrees -> radians -> rise-over-run
  
  i <- which(abs(slope) >= max_slope)
  
  conductance[i] <- 0
  
  # build matrix
  cm <- Matrix::spMatrix(nrow = n_cells, ncol = n_cells)
  
  cm[adj] <- conductance
  
  # include some spatial and other information to make
  # manipulating the network faster and easier
  # when passed to other functions
  # storing bounding box as vector to keep the size down
  bb8 <- terra::ext(x)
  bb8 <- c(bb8$xmin, bb8$xmax, bb8$ymin, bb8$ymax)
  
  xcrs <- terra::crs(x, describe = TRUE)
  xcrs <- paste0(xcrs$authority, ":", xcrs$code)
  
  tmin <- 1/max(conductance)
  tmax <- 1/min(conductance[conductance > 0])
  
  terrain <-
    list(
      "conductance" = cm,
      "neighbors"   = neighbors,
      "nrow"        = terra::nrow(x),
      "ncol"        = terra::ncol(x),
      "bb8"         = bb8,
      "crs"         = xcrs,
      "range"       = c("min" = tmin, "max" = tmax)
    )
  
  class(terrain) <- "terrain"
  
  return(terrain)
  
}

campbell <- function(x, decile = 50L) {
  
  if (length(decile) != 1) {
    
    stop("campbell accepts only one decile at a time.")
    
  }
  
  # campbell function wants degrees, not rise-over-run as in tobler
  # but, we're working with rr, so remember to convert radians to degrees with 180/pi!
  slope <- atan(x) * 180 / pi
  
  # values provided in Campbell 2019, Supplement, Table S1
  # for simplicity, I excluded all but the deciles
  deciles <-
    data.frame(
      decile = c(10L, 20L, 30L, 40L, 50L, 60L, 70L, 80L, 90L),
      a = c(-1.568, -1.71,  -1.858, -1.958, -2.171, -2.459, -2.823,  -3.371,  -3.06),
      b = c(13.328, 10.154,  8.412,  8.96,  10.064, 11.311, 12.784,  15.395,  16.653),
      c = c(38.892, 36.905, 39.994, 50.34,  63.66,  79.287, 98.697, 134.409, 138.875),
      d = c( 0.404,  0.557,  0.645,  0.649,  0.628,  0.599,  0.566,   0.443,   0.823),
      e = c(-0.003, -0.004, -0.004, -0.005, -0.005, -0.005, -0.005,  -0.005,  -0.0139)
    )
  
  ind <- which(deciles$decile == as.integer(decile))
  
  deciles <- deciles[ind, ]
  
  a <- deciles$a
  b <- deciles$b
  c <- deciles$c
  d <- deciles$d
  e <- deciles$e
  
  # lorentz distribution
  lorentz <- (1 / ((pi * b) * (1 + ((slope - a) / b)^2)))
  
  # modified lorentz
  (c * lorentz) + d + (e * slope)
  
}

#' Channel of reduced cost
#'
#' Reduces cost of travel along paths that cross through channels or
#' other aspects of a terrain that expedite travel. 
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#' @param channel an `sf` object.
#' @param .m a numeric between 0 and 1 applied as a weight to the cost of travel.
#'
#' @return a `terrain` representing cost of travel.
hf_channel <- function(x, channel, .m = 1) {
  
  rr <- terra::rast(
    nrow   = x$nrow,
    ncol   = x$ncol,
    extent = terra::ext(x$bb8),
    crs    = x$crs
  )
  
  channel <- terra::vect(channel)
  
  cells <- terra::cells(rr, channel)[, 2]
  
  adj <- cbind(
    "from" = as.integer(x$conductance@i + 1),
    "to"   = as.integer(x$conductance@j + 1)
  )
  
  # find incident edges, from or to cells intersecting barrier
  i <- which(adj[, 1] %in% cells | adj[, 2] %in% cells)
  adj <- adj[i, ]
  
  # want to reduce time, which means dividing conductance by .m
  x$conductance[adj] <- (x$conductance[adj] / .m)
  
  return(x)
  
}
primary_roads   <- roads |> filter(level == "P") |> st_buffer(dist = 30)
secondary_roads <- roads |> filter(level == "S") |> st_buffer(dist = 30)
tertiary_roads  <- roads |> filter(level == "T") |> st_buffer(dist = 30)

terrain <- dem |> 
  hf_terrain(max_slope = 35, decile = 30) |> 
  hf_channel(primary_roads,   .m = 0.2) |> 
  hf_channel(secondary_roads, .m = 0.5) |>
  hf_channel(tertiary_roads,  .m = 0.8)

remove(primary_roads, secondary_roads, tertiary_roads)
Code
#' Terrain to raster
#'
#' Converts a cost `terrain` into a raster representing
#' travel cost, specifically the average cost of moving from/to each
#' cell. Mostly useful for visualizing with `terra::plot()`
#' as a sanity check.
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#'
#' @return a `SpatRaster` with travel cost values.
hf_rasterize <- function(x) {

  rr <- terra::rast(
    nrow   = x$nrow,
    ncol   = x$ncol,
    extent = terra::ext(x$bb8),
    crs    = x$crs
  )

  # this part is basically what Jacob van Etten does with
  # gdistance::raster(TransitionLayer), though that offers more options
  # here I am just averaging the vertical and horizontal sums
  total_row <- Matrix::rowSums(x$conductance)
  total_col <- Matrix::colSums(x$conductance)

  #logical sparse matrix, identifies number of adjacent cells
  lm <- methods::as(x$conductance, "lMatrix")

  n_row <- Matrix::rowSums(lm)
  n_col <- Matrix::colSums(lm)

  # average cost of moving from/to each cell
  # values <- total_cost / n_adjacent
  v_row <- total_row / n_row
  v_col <- total_col / n_col

  values <- (v_row + v_col)/2

  terra::setValues(rr, (1/values))

}

gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hf_rasterize(terrain), xy = TRUE),
    aes(x, y, fill = lyr.1)
  ) +
  scale_fill_viridis(
    name = "Seconds",
    breaks = c(10, 30, 50)
  ) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y"), 
    expand = expansion(add = 1000)
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = 3.8, 
      barheight = 0.6,
      raster = FALSE,
      frame.linewidth = 0.2,
      ticks.linewidth = 0.7,
      title.position = "top"
    )
  ) +
  theme(
    legend.position = c(0.03, 0.98),
    legend.direction = "horizontal",
    legend.spacing.x = grid::unit(0.05, "cm"),
    legend.text = element_text(size = rel(0.5), margin = margin(l = 2)),
    legend.title = element_text(
      size = rel(0.8), 
      margin = margin(b = 2, t = 2),
      hjust = 0
    )
  )

ggsave(
  plot = gg,
  filename = here("figures", "cost-surface.png"),
  dpi = 600,
  width = 3.5,
  height = 3.5 * (dy/dx)
)

4.1 To springs

Function to calculate paths and travel times is in the folded code chunk below.

Code
#' Shortest path between points
#'
#' Calculate the "shortest" or least-cost paths between multiple start and end points.
#' For each from-point, the shortest path to every to-point is calculated.
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#' @param from an `sf` specifying start locations, must be POINT or MULTIPOINT.
#' @param to an `sf` specifying end locations, must be POINT or MULTIPOINT.
#' @param add_cost logical, whether to include a column of travel costs.
#'
#' @return An `sf` object with least-cost `LINE` features and two columns: from and to.
#'   Values in the from- and to-columns are `1:nrow(from)` and `1:nrow(to)`, respectively.
#'   If `add_cost = TRUE`, the `sf` object also includes a travel cost column.
hf_hike <- function(x, from, to, add_cost = FALSE) {
  
  from_xy <- sf::st_coordinates(from)[, 1:2, drop = FALSE]
  to_xy   <- sf::st_coordinates(to)[, 1:2, drop = FALSE]
  
  rr <- terra::rast(
    nrow   = x$nrow,
    ncol   = x$ncol,
    extent = terra::ext(x$bb8),
    crs    = x$crs
  )
  
  from_cells <- terra::cellFromXY(rr, from_xy)
  to_cells   <- terra::cellFromXY(rr, to_xy)
  
  graph <- igraph::graph_from_adjacency_matrix(
    x$conductance,
    mode = "directed",
    weighted = TRUE
  )
  
  # invert conductance to get travel cost
  igraph::E(graph)$weight <- (1 / igraph::E(graph)$weight)
  
  # return a list of lists of vectors of vertices on paths
  shorties <-
    lapply(
      from_cells,
      function(z) {
        
        res <-
          igraph::shortest_paths(
            graph,
            from = z,
            to = to_cells,
            mode = "out"
          )
        
        res$vpath
        
      })
  
  # collapse (list of lists of vectors) to (list of vectors)
  shorties <- unlist(shorties, recursive = FALSE)
  
  # now we can just make a bunch of linestrings out of them
  line_list <-
    lapply(
      shorties,
      function(z) {
        
        cells <- as.vector(z)
        
        # to handle case where from-vertex == to-vertex
        if (length(cells) == 1) cells <- rep(cells, 2)
        
        xy <- terra::xyFromCell(rr, cells)
        
        sf::st_linestring(xy)
        
      }
    )
  
  # and then an sf
  sf_col <- sf::st_sfc(line_list, crs = x$crs)
  
  short_paths <- sf::st_sf(geometry = sf_col)
  
  # add indices
  short_paths <- transform(
    short_paths,
    from = rep(1:nrow(from_xy), each = nrow(to_xy)),
    to   = rep(1:nrow(to_xy), times = nrow(from_xy))
  )
  
  if (add_cost) {
    
    cost <- igraph::distances(graph, from_cells, to_cells, mode = "out")
    
    short_paths <- transform(short_paths, cost = c(t(cost)))
    
  }
  
  return(short_paths)
  
}
to_springs <- hf_hike(
  terrain, 
  from = clusters, 
  to = springs,
  add_cost = TRUE
)

to_springs <-
  to_springs |>
  group_by(from) |> 
  slice_min(cost, n = 2) |> # for each cluster, get two lowest cost paths
  mutate(cost = (cost/3600)) |> # convert to hours 
  summarize(
    cluster = unique(from),
    spring = paste(to, collapse = ","),
    cost = mean(cost)
  ) |> 
  select(cluster, spring, cost)

4.2 To cities

Please note that travel time is not strictly to cities, but to the edge of the project area. As noted above, for the northern watershed, the time to La Paz equals the time it takes to leave the watershed plus half an hour.

north <- clusters |> 
  st_filter(
     st_union(filter(watersheds, id %in% c(3294, 3296, 3299)))
  )

south <- clusters |> 
  st_filter(
    filter(watersheds, id == 3326)
  )

### NORTH ###
out_north <- hf_hike(
    terrain, 
    from = north, 
    to = terminus[1, ], 
    add_cost = TRUE
  ) |> 
  mutate(
    cluster = north$cluster,
    city = "Ciudad Constitucion", 
    cost = (cost/3600),
    watershed = "north"
  )

### SOUTH ###
out_south <- hf_hike(
    terrain, 
    from = south, 
    to = terminus[2, ], 
    add_cost = TRUE
  ) |> 
  mutate(
    cluster = south$cluster,
    city = "La Paz",
    cost = (cost/3600),
    watershed = "south"
  )

to_cities <- 
  bind_rows(out_north, out_south) |> 
  mutate(time_to_paz = ifelse(watershed == "north", cost + 0.5, cost)) |>
  select(cluster, city, cost, time_to_paz)

remove(north, south, out_north, out_south)

4.3 LCP Map

Code
od_points <- bind_rows(
  od_points |> mutate(destination = ifelse(variable == "spring", "To springs", "To cities")),
  od_points |> filter(variable == "cluster") |> mutate(destination = "To springs")
)

paths <- bind_rows(
  to_springs |> mutate(destination = "To springs"),
  to_cities |> mutate(destination = "To cities")
)

gg <- ggplot() +
  geom_raster(
    data = as.data.frame(hillshade, xy = TRUE),
    aes(x, y, fill = shade)
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide = "none"
  ) +
  new_scale_fill() +
  geom_raster(
    data = as.data.frame(dem, xy = TRUE),
    aes(x, y, fill = elevation),
    alpha = 0.7
  ) +
  scale_fill_distiller(
    palette = "Greys", 
    na.value = "transparent",
    guide="none"
  ) +
  new_scale_fill() +
  geom_sf(
    data = watersheds, 
    fill = "transparent", 
    color = alpha("darkblue", 0.35),
    linewidth = 0.2
  ) +
  geom_sf(
    data = roads, 
    color = "black",
    linewidth = 0.3
  ) +
  geom_sf(
    data = paths |> st_buffer(125) |> group_by(destination) |> summarize(),
    fill = "#C7FFED",
    color = "#00291D",
    linewidth = 0.2
  ) +
  geom_sf(
    data = od_points, 
    aes(shape = variable, fill = variable, color = variable),
    stroke = 0.3,
    size = ifelse(od_points$variable == "city", 2.2, 1.5)
  ) +
  facet_wrap(vars(destination)) +
  scale_color_manual(values = c("grey15", "white", "white")) +
  scale_fill_manual(values = alpha(c(point_colors, 0.75))) +
  scale_shape_manual(values = c(23, 22, 21)) +
  coord_sf(datum = 4485) +
  scale_x_continuous(
    breaks = pretty_breaks(bb8, "x"), 
    expand = expansion(add = 1000)
  ) +
  scale_y_continuous(
    breaks = pretty_breaks(bb8, "y", n = 3), 
    expand = expansion(add = 1000)
  ) +
  theme(
    legend.position = c(0.01, 0.98),
    strip.background = element_blank(),
    strip.text = element_text(size = rel(1.1), hjust = 0)
  ) 

ggsave(
  plot = gg,
  filename = here("figures", "least-cost-paths.png"),
  dpi = 600, 
  width = 5.75,
  height = 5.75 * (dy/(2*dx)) + 0.47
)

4.4 Save

write_sf(
  to_springs,
  dsn = gpkg,
  layer = "path_to_springs"
)

write_sf(
  to_cities,
  dsn = gpkg,
  layer = "path_to_cities"
)

5 References

Campbell, Michael J., Philip E. Dennison, Bret W. Butler, and Wesley G. Page. 2019. “Using Crowdsourced Fitness Tracker Data to Model the Relationship Between Slope and Travel Rates.” Applied Geography 106: 93–107. https://doi.org/https://doi.org/10.1016/j.apgeog.2019.03.008.

6 Session Info

Code
# save the session info as an object
pkg_sesh <- sessioninfo::session_info(pkgs = "attached")

# inject the quarto info
pkg_sesh$platform$quarto <- paste(
  quarto::quarto_version(), 
  "@", 
  quarto::quarto_path()
  )

# print it out
pkg_sesh
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16 ucrt)
 os       Windows 11 x64 (build 22621)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  English_United States.utf8
 ctype    English_United States.utf8
 tz       America/Denver
 date     2024-01-23
 pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
 quarto   1.4.545 @ C:\\PROGRA~1\\Quarto\\bin\\quarto.exe

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 dplyr       * 1.1.4   2023-11-17 [1] CRAN (R 4.3.1)
 ggnewscale  * 0.4.9   2023-05-25 [1] CRAN (R 4.3.1)
 ggplot2     * 3.4.4   2023-10-12 [1] CRAN (R 4.3.2)
 here        * 1.0.1   2020-12-13 [1] CRAN (R 4.3.1)
 igraph      * 1.6.0   2023-12-11 [1] CRAN (R 4.3.2)
 sf          * 1.0-15  2023-12-18 [1] CRAN (R 4.3.2)
 terra       * 1.7-65  2023-12-15 [1] CRAN (R 4.3.1)
 viridis     * 0.6.4   2023-07-22 [1] CRAN (R 4.3.1)
 viridisLite * 0.4.2   2023-05-02 [1] CRAN (R 4.3.1)

 [1] C:/Users/kenne/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.1/library

──────────────────────────────────────────────────────────────────────────────